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A many-flavor electron gas (MFEG) in a semiconductor with a valley degeneracy ranging between 
6 and 24 was analyzed using diffusion Monte Carlo (DMC) calculations. The DMC results compare 
well with an analytic expression derived by one of us [Phys. Rev. B 78, 035111 (2008)] for the 
total energy to within ±1 % over an order of magnitude range of density, which increases with 
valley degeneracy. For BigTea (six-fold valley degeneracy) the applicable charge carrier densities are 
between 7 x 10^® cm~^ and 2 x lO^'^cm"^. DMC calculations distinguished between an exact and 
fT^ ' a useful approximate expression for the 24-foId degenerate MFEG polarizability for wave numbers 

2pp < g < 7pp. The analytical result for the MFEG is generalized to inhomogeneous systems by 
means of a gradient correction, the validity range of this approach is obtained. Employed within 
' ^ ' a density-functional theory calculation this approximation compares well with DMC results for a 

,1 ^ I quantum dot. 

» . 

. PACS numbers: 71.15.Mb,71.10.Ca,02.70.Ss 

B 

^ • I. INTRODUCTION 

o . 

^ _ Good quantum numbers, that describe conserved quantities as a quantum system evolves, derive their significance 
from their connection to the powerful conservation laws of physics. In addition to the famiUar examples of spin 
and crystal momentum, under some circumstances electrons in solids can have an additional quantum number that 
" ^ • distinguishes them, which we call the flavor; we denote the total number of flavors by v. One example of such a system 
■ are semiconductors and scmimetals that have degenerate conduction-band valleys, the flavor denotes the electron's 
, valley. Examples of multi- valley semiconductors include Ge, which as shown in Fig. [1] has four degenerate valleys 

• (N.B. not eight, as valleys at the Brillouin zone vertices overlap). Si has six degenerate valleys, a Ge-Si alloy has ten 
] degenerate valleys, and Pbi-x-j/Sn^jMnyTe has twelve valleys in the S band The system has been experimentally 

• realized as an electron- hole liquid that forms in drops 0, Q . In these systems the number of flavors (the number of 
T— I ■ valleys) is well defined and there are strong Coulomb interactions between particles which motivates the analysis. This 
OO . is in contrast to several other systems in which the number of flavors is poorly defined such as heavy fermions [1, H, @1 j 

' charged domain walls 0] , a super-strong magnetic field 0] , and spin instabilities [1, ; or where the number of flavors 
J> . is well deflned but interactions between particles are weak such as ultracold atoms in optical lattices [ll|, [l^, ■ 

The properties of a many-flavor electron gas (MFEG) in a semiconductor were first studied analytically for the 
normal phase by Andryushin ct al. 0, and for the superconducting phase by Cohen 14 1. Recently one of us Q 
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■ extended the MFEG analysis by finding an energy functional and gradient expansion, which allowed the study of 
inhomogeneous systems. However, the analytical treatment was limited to consider the same contributions to the 
energy as in the random phase approximation (these contributions dominate in the many-flavor limit). To go further 
requires numerical calculations, the only example of which for a MFEG to date lla used a self-consistent approach 
for the local fleld correction formulated by Singwi et al. [31 (STLS), see also Ref. jl7l). The method was later appHed 
to charge impurities by Bulutay et al. The calculations of Ref. [lH| were performed for v <(S, too few flavors to 
gauge the applicability of the analytic many-flavor approximation, which is estimated to apply at around six or more 
flavors d]. 

In this paper we follow the suggestion of Gold [l^ , and present the results of what are expected to be more accurate 
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FIG. 1: The Ge band-structure in the [111] direction calculated using a plane-wave pseudopotential method [T^. The Fermi 
energy is at i5 = OeV; below are valence bands with the holes centered around H, above are conduction bands. The first 
conduction band valley is highlighted in bold, low-lying conduction-band electrons are centered around C. 



diffusion Monte Carlo (DMC) [13, [HI, [H, [H] calculations on the MFEG for v < 24, which should allow us to verify 
the analytical MFEG approach. We then examine aspects of the many-flavor approximation that have not yet been 
studied computationally: in Sec. IIVI we compare the analytical density-density response function derived in Sec. II Bl 
with that predicted using DMC. Once verified this allows us in Sec.|V]to employ a gradient expansion within density- 
functional theory (DPT) to find the ground state of a quantum dot, we compare results with DMC calculations and 
examine the validity of the gradient expansion. 

We adopt the atomic system of units: that is e'^ ~ h ~ m ~ l/(47reo) = 1. The mass m — mem* is defined to be 
the electron mass, TOo, multiplied by a dimensionless effective mass m* appropriate for the conduction-band valleys, 
which when m* = 1 will recover standard atomic units. We assume the valleys all have the same dispersion profile 
and so the same effective mass, Andryushin et al. [H outlined a method of calculating a scalar effective mass for 
anisotropic valleys. With the above definitions, energy is given in terms of an exciton = E\jn* ^ where i^h is the 
Hartree energy, and length = a^/m* in terms of the Bohr radius cq. To denote density we use both the number 
density of conduction-band electrons n and the Wigner-Seitz radius r^. 

Before presenting the numerical results, to orient the discussion, we describe the basic physics of the MFEG and 
review the analytical results of Ref. [3] that will be computationally verified in this paper. 



A. Introduction to a MFEG 



In a low temperature MFEG, the number of fiavors v, number density of conduction-band electrons n, and Fermi 
momentum p-p are related through 

At fixed electron density, the Fermi momentum reduces with increasing number of flavors as pp oc v~^/'^ ^ so each 
Fermi surface encloses fewer states. The semiconductor hole band-structure often has a single valence-band minimum 
at the F point, such as in Ge, see Fig.[l] hence we assume the holes are heavy and are uniformly distributed, providing 
a jellium background. 

For a constant number density of particles, the density of states at the Fermi surface, g, rises with increasing number 
of flavors as g oc v^/Ep oc v"^^^. Therefore, the screening length estimated with the Thomas-Fermi approximation 
is = (47re^g)^^/^ oc i'^^^^, and the ratio of the screening to Fermi momentum length-scale varies with number 
of flavors as pp/n oc v'"^!^ . In the many-fiavor limit v \^ the screening length is much smaller than the inverse 
Fermi momentum, kT^ <C Pp^, and so the dominant electron-electron interactions have characteristic wave vectors 
which obey q pp. This is in direct contrast to the random phase approximation (RPA) where pp 3> n, although 
in both the many-flavor and the RPA, the same Green function contributions with empty electron loops dominate 
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diagrammatically [2, These diagrams contain the greatest number of different flavors of electrons, and as 3> 1 
therefore have the largest matrix element. Since q ^ Pf, the typical length-scales of the MFEG are short, this 
indicates that a local density approximation (LDA) could be applied. This motivation is in addition to the usual 
reasons for the success of the LDA in DFT [251 , namely that the LDA exchange-correlation hole need only provide a 
good approximation for the spherical average of the exchange-correlation hole and obey the sum rule [1^ . 



B. Polar izability 

In the many-flavor limit the exact result for the polarizability of a MFEG at wave vector q, and Matsubara frequency 
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This quantity governs the density-density response of the MFEG so is important to verify. Since Eqn. ([3]) has a simple 
form it can be used to calculate further properties of the MFEG such as homogeneous energy in Sec. II Gl and the 
gradient expansion in Sec. II Di which further motivates its numerical verification. 



C. Homogeneous energy 

Starting from the approximate expression for polarizability, Eqn. ([3]), it can be shown that the total energy of a 
MFEG, including all the exchange and correlation contributions is Q 

^'int 

where Ago = r(-5/4)r(3/4)/(27r^/4)(£;*a*^/'') w 0.574447(£;i*a*^/*) and Ei^^ denotes the interacting energy (which 
would be zero if electron-electron interactions were ignored) . 

In Ref. 3 it was suggested that this relation for the total energy applies over a density range, at 99 % accuracy, 
Q.Oiv <^ nag'^ <^ (0.0741^)"*, which widens with number of flavors as (see also Ref. [2[). Gonsidering the number 
of flavors where the range of validity vanishes indicates that the many-flavor limit will apply if there are ten or more 
flavors. An alternative estimate for the density range is found in Sec. IIV Al bv comparing the analytical result with 
DMG calculations. 



D. Gradient correction 



The applicability of the LDA in a MFEG motivates the search for a gradient expansion to the energy Eqn. (|4]) 
as a way to analyze inhomogeneous systems such as electron-hole drops and quantum dots. The typical momentum 
transfer in the MFEG is q ^ 4(fi,aQ ^^^)ri^/^, which defines the shortest length-scale over which a LDA can be made, 
therefore, the maximum permissible gradient in electron density is |Vn|niax ^ ~ 4(na*"^/*)n5/4. A gradient 
expansion will break down for phenomena with short length-scales, for example mass enhancement [28j . If electron 
density is smoothly varying then starting from Eqn. the gradient correction to the energy for a MFEG is Q 

E.E^.i'y^, (5) 

8 n 

where Eq is the energy of a homogeneous MFEG with density n, see Eqn. As discussed in Sec. II Al this gradient 
expansion would be useful for DFT calculations and so its computational verification is important. 
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II. COMPUTATIONAL METHOD 

111 this section we briefly describe the two computational methods that we used, variational Monte Carlo (VMC) 
and Diffusion Monte Carlo (DMC) [l^. These are quantum Monte Carlo (QMC) methods, chosen since DMC gives 
the exact ground state energy subject to the fixed node approximation, and both are expected to give more accurate 
results than the STLS approach used by Gold [ll] . 

The VMC method uses a normalizable and difFcrentiable trial wave function of the form discussed below. The 
Metropolis algorithm p9l | is used to sample the wave function probability density l^'xP using a random walk, and make 
an estimate of the local energy £'L(r) = ^'T(r)~"'^-ff4'T(r). In order to obtain the ground state one could minimize 
the spatial average of the local energy with respect to the free parameters in the trial wave function. However, it 
is computationally more stable to minimize the variance in the estimates of the local energy. As VMC obeys the 
variational principle by construction, it yields an upper bound to the true ground state energy. 

The more accurate DMC algorithm is a stochastic method that begins with a trial or guiding wave function, in this 
case the optimized VMC trial wave function. The DMC method is based on imaginary time evolution, which when 
using the operator e"*^^^^"^-* projects out the ground state wave function from the trial wave function, and yields an 
estimate of the ground state energy, E^. The nodal surface on which the wave function is zero (and across which it 
changes sign) is fixed (2ll.[30l| to be that of the trial wave function, this ensures that the fermionic exchange symmetry 
is maintained. The DMC algorithm produces the exact ground state energy subject to the fixed node approximation, 
and is also variational so gives an accurate upper bound to the true ground state energy once the population control 
bias and finite time-step bias are eliminated. The algorithm used closely follows that described in Rcf. (31i |. 

In our QMC calculations we use a Slater-Jastrow [l^, [s^, trial wave function. The Slater part of the wave 
function is a product of determinants, each one corresponding to a different electron spin or flavor. Each determinant 
is over the spatial orbitals of electrons occupying the lowest energy levels. The determinant changes sign when rows or 
columns are swapped, this ensures that the wave function is antisymmetric under exchange of electrons with the same 
flavor and spin. The Slater wave function itself is not the ground state of an interacting electron gas, so to improve 
the wave function, variational degrees of freedom that account for two-body correlations are included within a Jastrow 
factor. The Jastrow factor is symmetric under particle exchange so does not alter the particle exchange symmetry 
of the wave function. Furthermore, the Jastrow factor is always positive so docs not alter the wave function nodal 
surface. The Jastrow factor contains a two-body polynomial term uijij) = F{rij) X]f=2 '^'''ij' ^ power series form (ssj 
in electron separation Vij with optimizable parameters, ai. The term F{rij) ensures that the Kato cusp conditions are 
satisfied [13]. To ensure that electron-electron correlations do not extend beyond the simulation cell, the term is cutoff 
at the Wigner-Seitz radius. To treat longer-ranged correlations, the Jastrow factor includes a two-body plane-wave 
expansion, p(ry) = cos(Ga ■ Those reciprocal lattice vectors, {G^i}, that are related by the point 

group symmetry (denoted by A) of the Bravais lattice share the same optimizable parameters, a a- To ensure accuracy 
we checked the stability of the VMC results when the expansion order of the u and p terms was increased. At all 
densities the Jastrow factor optimized cutoff lengths took the maximum allowed value (the Wigner-Seitz radius). 

The DMC calculations were performed with 57 different reciprocal lattice vectors and, following Ortiz and Ballone 
[ssf . Ceperley [s^, and Ceperley and Alder [s^], further VMC calculations were performed at other system sizes 
(27, 33, 57, and 81 reciprocal lattice vectors) to derive the parameters to extrapolate the DMC energy to infinite 
system size. Additionally, all the DMC results were extrapolated to have zero time-step between successive steps in 
the electron random walk. In DMC simulations the acceptance probability of a proposed step in the random walk 
exceeded 99%. We used 300 DMC configurations, comparable to the 200-300 used by Ortiz and Ballone [s^], and 
checked for population control bias by ensuring that ground state energy estimates did not vary with a changing 
number of configurations. All the QMC calculations were performed using the CASINO computer program [s^. 

III. HOMOGENEOUS MFEG 

We start with the simplest possible system to analyze numerically, the homogeneous MFEG, this provides not only 
a suitable system to validate both theory (Sec. II C|) and the QMC many-flavor calculations, but should also confirm 
the range of densities over which the many-flavor approximation applies. The 3D homogeneous electron gas [v = 1) 
has been studied before using QMC [1^, HI, [131 and these studies provide a useful guide to the method we should 
follow. 

To calculate the interaction energy E^at we subtracted the theoretical Thomas-Fermi kinetic energy from the DMC 
ground state energy (sec Eqn. ([5])). At each of 6, 12, 18, and 24 flavors we performed flve DMC calculations and 
interpolated to flnd where theory and DMC results agree to within ±1%. Results in Fig. [2] show that for v > Q the 
theory applies over at least an order of magnitude in density to an accuracy of ±1 % - the theory can be applied at 
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FIG. 2: The lower panel shows the fractional difference of DMC interaction energy Eumc from the model Ei-at with MFEG 
density n (and Wigner-Seitz radius r^) for different numbers of flavors, the dotted lines show ±1% disagreement. The central 
panel bars highlight the numerical region of applicability, the gray shaded area represents the analytically predicted region of 
±1% applicability. The upper panel shows the total energy for 6, 12, 18 and 24 flavor electron gases. 



fewer flavors than expected. For fewer than ~ 12 flavors the valid logarithmic range of the theory increases with z/, 
the 18 and 24 flavor resuhs show a dramatic increase in the range of validity, especially on the high density side. In 
the limit of many flavors {v > 12) the expected 99% range of validity O.OSz^ <C na^'^ <C (0.074j/)'* is approximately 
consistent with the computationally predicted ±1 % region, therefore the minimum number of flavors required for aU 
aspects of the many-flavor theory to be valid is approximately ten. 

For Si with m* = 1.08 the many-flavor limit applies to an accuracy of ±1% for a charge carrier concentration 
between 4 x 10^'^ cm^'^ and 1 x 10^^ cm^'^, this is greater than the typical maximum carrier density ^ 1 x 10^^ cm^"^ 
and so in Si the formalism is not applicable. In systems with a low effective mass, for example the ly = 6 material 
Bi2Te3 used in thermoelectric cooling, which has m* = 0.06 [1^, SO, |4l| , the required charge carrier concentration 
is between 7 x lO^^cm"^ and 2 x 10^°cm~^, which compares favorably with the typical maximum carrier density 
~ 1 X 10^^ cm~'^ and so the many-flavor limit formalism could be applied to low effective mass materials. 

The STLS results of Gold [1^1 sA, v = 1 were ^ 3.4% less negative than the DMC results of Ortiz and Ballone [s^, 
and at J/ = 6 were ~ 3.1% less negative than our DMC results. This represents a significant difference between our 
and the STLS results when looking for the 1% range of validity, highlighting the need for the more accurate DMC 
calculations. The range of validity at ±1 % up to at least 24 flavors is to the high density side of the minimum in 
the total energy seen in Fig. [2l but the minimum nmin oc u^^^ lies within the region of validity for higher v. 7i,nin 
is the density expected to be seen in physical systems such as electron- hole drops, the good agreement of the theory 
with DMC results at this density indicates that the theory could be usefully applied to investigate the properties of 
physical systems, see for example Ref.Q. 
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IV. STATIC DENSITY-DENSITY RESPONSE 



Having verified the homogeneous system behavior we may now proceed and computationaUy examine inhomo- 
geneous behavior through the static density-density (hnear) response function Eqn. ([3]). The polarizability is an 
important quantity used 0] to develop both homogeneous theory and the gradient correction, the density-density re- 
sponse function itself also governs the electrical response properties, for example polarization, screening, and behavior 
in an external potential; it is therefore useful to verify this response before applying the theory to model systems. We 
examine l/e(q), the quantity probed experimentally (4^ . 

DMC has previously been used to find the static density-density response of single-fiavor systems: Sugiyama et al. 
[1^ applied the method to charged bosons, the density-density response of the electron gas was calculated by Moroni 
et al. [31 (in two dimensions), and Bowen et al. [1^1 and Moroni et al. [iHl (three dimensions). However, density- 
density response has not been studied numerically in a many-flavor system. Here we employ two methods to find 
the density-density response function. The more accurate and computationally efficient method of calculating the 
response is to examine the ground state energy, calculated using DMC. A VMC energy based estimate and an estimate 
using the induced electron density are used to check the accuracy of the trial wave function. 

Before the results are described in Sec. IIVCI we outline the theory behind the two methods used to estimate 
the response, firstly in Sec. IIV Al bv using the ground state energy variation, and secondly in Sec. IIV Bl through the 
magnitude of the periodic density modulation. 



A. Ground state energy variation 



To calculate the density-density response we use a weak probe so that the density response is solely due to the 
properties of the homogeneous system. We apply a static {uj = 0) monochromatic perturbative external potential 
U{y) = ?7qC0s(q • r) to the homogeneous MFEG, corresponding to the background charge having an additional 
sinusoidal variation noxt(r) = 7iqCos(q • r). The external potential and external charge are linked [i^ through 
Poisson's equation by 

"oxt(k) = %^((5k.q + <5k.-q) . (6) 

We assume that different Fourier components are independent, the density response to an external potential with 
wave vector q and frequency m is only at that wave vector and frequency so the induced charge is nind(k) = 

{{^i^lUq ~ (^k)o) (i^k.q + (^k.-q). Here {h-k)uci expectation value of the charge density Fourier component 

at wave vector k with an applied external potential C/q, and (?T.k)o is the same but in the homogeneous case with 
no external potential. Linear response theory gives the static density-density response function as the ratio of the 
induced charge density and the perturbing external charge density so 

^ = 1 + ^ (K>.. - . (7) 

If the external potential is small relative to other typical energies the density response is determined solely by the 
properties of the homogeneous MFEG. We can expand in small so that 



(nk)f/„ - ("k)n ~ U, 



d(nk) 



q 



dt/o 



(8) 

Ua=0 



where the induced charge density is calculated by considering the dependence of the ground state energy E on the 
magnitude of the external field. Substituting this into Eqn. gives an expression for the density-density response 



1 Stt d^E 

= 1 + ^ 



e(q) g2 d[/2 



(9) 

(7a =0 



To recover the density-density response function at a particular wave vector, several QMC calculations were per- 
formed at that wave vector for different amplitudes of the external field. A polynomial fit was made to the ground 
state energy so as to extract the second derivative. To investigate the lowest order polarizability the applied external 
field should be as small as possible yet still give statistically significant results, to ensure this we checked that the 
ground state energy showed only quadratic behavior with applied field amplitude. A further convenient way to check 
the perturbing field is sufficiently small is to ensure the electric field of the external potential is less than the typical 
electric field strength between two neighboring electrons, e/r^. 
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B. Induced charge density measurement 

As the external potential is perturbative we use the same plane-wave basis set as employed for the calculations on 
the homogeneous MFEG described in Sec. IIIII To account for the modulating density, following Moroni ot al. [i^ . 
Bowen et al. [45|, and Moroni ct al. [i^ we introduce a new q term into the Jastrow factor of the form 

(7(r,) = fecos(q • Ti) , (10) 

where b is an optimizable parameter, the position of the ith electron, and the wave vector q corresponds to that of 
the perturbative external potential. As b is small, the charge density induced by the perturbative external potential 
is nind ~ 26cos(q • r^). From Eqn. ([6]) and Eqn. ([7|) it follows that 

1 87r6 

1 + -^- (11) 



e(q) q'Uq 

The optimized value of b was found by variance minimization during a VMC calculation. The relationship then allows 
us to derive an estimate for the density-density response function for each separate C/q, typically four values were 
averaged to give a final estimate for the density-density response. 



C. Results 



We chose to find the polarizability for a MFEG with v = 24 and Tg = 0.6ag. This lies at the lower bound of the 
range of validity near to the minimum in the energy (see Fig. [2]) at a density expected to be seen in physical systems. 
This density was also chosen since it had most of the polarizability curve 0.25 < 1/e < 1 in the region of applicability 
q > 2pp. Boundary conditions mean that the external potential must be periodic over the simulation cell, therefore 
the external potential wave vector q must be a reciprocal lattice vector. We checked that if the Jastrow factor q term 
wave vector was changed so that it was incommensurate with the external potential then following optimization & = 
within statistical errors; this verified the linear response assumption that Fourier components arc independent. 

The results of the calculation are shown in Fig. [31 The DMC results obtained by considering the variation in ground 
state energy (see Sec. IIV Ap better fits the exact than approximate expression for the polarizability, and though error 
bars are large can distinguish between the two within one standard deviation. This shows that QMC results can exceed 
the accuracy of the approximation made in Eqn. ([3]) , though that estimate remains useful. The positive agreement 
verifies the theory and confirms the accuracy of the CASINO simulations. 

The ground state energies calculated by VMC were used in the same way as the DMC results to find the density 
response and provide a reasonable fit, though here error bars are large so comparison is difficult. Following the 
prescription in Sec. IIVBI we also derived values for the density-density response function using the charge density 
modulation at the wavelength of the perturbing potential, C/q. These values agreed within statistics though carried 
a larger uncertainty than those derived using the ground state energy. Both of these alternative methods appear to 
underestimate the density-density response. These results are consistent, a smaller charge density response gives a 
smaller coefficient in the Jastrow factor q term and a smaller reduction in ground state energy. Nevertheless, the 
reasonable agreement of both VMC estimates and to the DMC results indicates that the trial wave function had an 
adequate nodal surface. 



V. GRADIENT CORRECTION 



It was important to verify the density-density response as it is a key component to the many-flavor formalism and 
could be applied to other many- flavor systems where density is expected to be inhomogeneous, for example junctions 
and the response to defects and impurities. Now that it has been verified, we may proceed to consider a quantity 
derived from it: the gradient expansion, Eqn. ([5]), which is also useful for analyzing systems with inhomogeneous 
density. Once we have investigated the validity of such an expansion we can apply the formalism to quantum dots, 
chosen since they have a large controllable variation in electron density so should provide a good test of the gradient 
expansion. Quantum dots are commonly made in many-flavor semiconductor materials so can be modeled using a 
many-flavor formalisrn, and are a system in which there is current research interest. 

Quantum dots [ZtI . [isj have not previously been studied in the many-flavor limit though there have been several 
previous computational studies of a single- flavor electron gas confined in a quantum dot. Previous QMC simulations of 
quantum dots include Pollock and Koch [i^, Harju ct al. [l^l performed VMC calculations for parabolically confined 
electrons in circular dots. Bolton [5l| performed fixed-phase DMC simulations. Path-integral QMC calculations have 
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FIG. 3: The density-density response 1/e versus the wave vector |q| of a MFEG with i/ = 24 and = O.Gag. The sohd curve 
shows the exact resuh l/texact (Eqn. (O), the dotted curve the Eqn. (j3} approximation. The shaded grey region |q|/2 < pp is 
where the many-flavor limit breaks down. The points show the values for the permittivity calculated from QMC results, the 
circle is from modulated charge measurements (CM), the triangle from VMC energy and the cross from DMC energy. The 
lower panel plots the actual response, the upper panel shows the deviation of response from the exact theoretical result Eqn. ([2]) 
with standard error bars. 

also been performed [s^. [ssl . [53 . these showed poor agreement with results from exact diagonalization [5^1. Benedict 
et al. [IHl, Williamson et al. |57l|. Puzdcr ct al. [HHl all compared the optical band-gap between DMC calculations 
and results from other methods. For circular quantum dots Pcdcriva ct al. [s^ found the ground-state using both 
DMC, a local spin density approximation method, and Hartree-Fock, they then directly compared the ground-state 
energy, correlation energy, and spin density profiles. Ghosal et al. [iOl also used DMC to investigate circular quantum 
dots. Quantum dots have successfully been investigated using DFT [H, [6l|, [63, [gl], Pederiva et al. [5^ found the 
local spin density approximation method predicted ground-state energies that were typically 2 % greater than DMC 
energies, Ferconi and Vignale [gH obtained a 3 % agreement between currcnt-density-functional theory and exact 
diagonalization results. 

A. Method 

Before describing the study of quantum dots using a many-flavor functional in detail we first outline the general 
strategy of the numerical calculations. Firstly, a DFT calculation using the many-flavor functional (including the 
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gradient approximation) was performed using a plane-wave basis set. This produced an estimate of the ground-state 
energy and density according to the many-flavor theory. It also provided a trial wave function that was converted 
to a B-spline basis set and, with Jastrow factor, was optimized in a VMC calculation, in preparation for a DMC 
calculation. Finally, the DMC calculation gave a second estimate of the ground state energy and density, exact only 
for the fixed node approximation. This estimate was compared with the DFT calculation, and also gave an insight 
into the accuracy of the many-flavor theory. 

Here we carried out simulations on a quantum dot with a harmonic external potential of the form V = where 
r is the distance to the center of the quantum dot containing a MPEG with 12 flavors. This potential was chosen as it is 
simple, continuous, reahstic [13, [11], and has been used in previous computational studies [sol. [sol. [6ll. [6^. [63l. [66l. [btI. [68| . 
Filled shells in this potential correspond to 1, 4, 10, 20, 35, . . . orbitals (whose degeneracy may be reduced by electron- 
electron interactions). In DFT we used a supercell containing a single dot to model the aperiodic system with periodic 
boundary conditions, in DMC non-periodic calculations with just a single quantum dot were performed. The cubic 
cell was large enough that the trial wave functions had reduced by at least a factor of 10~^ at its boundary. 

Trial wave functions were generated using the DFT program 3Ddotdft, an extended version of DOTDFT [g^. This 
used the many-flavor functional with gradient approximation so had energy density 

eKr)) = -A3Dn5/4 + e^. (12) 

8n 

A new parameter ^ was introduced that multiplies the gradient term, which allowed us to adjust its size; C = 1 gives 
the correct analytical expression, and ^ = the functional without a gradient expansion. 

The VMC simulations, run in CASINO, used a B-spline basis set [tO, [zH because a localized basis set offers 
significant performance advantages over plane-waves. The wave function was optimized in VMC with a Jastrow 
factor containing the two-body polynomial u term and two-body plane-wave term p with the same form as used in 
Sec, mil and Sec. llVi and a one-body electron-potential term x(''i) = Sm=2 /^m'^i" with F determining behavior 

at the cutoff length, ri the distance of the ith electron from the center of the potential, and the (3m being optimizable 
parameters; we also note that the x term has no central cusp. 

The many-flavor functional incorrectly adds in the self-interaction energy of each electron to its own Coulomb 
potential. One way to correct for this is to add an additional term to the density functional [zl, [ll]. However, as the 
number of flavors is increased the ratio of the correct interaction (oc v^) to incorrect self-interaction (oc v) increases 
as ~ 2z/ — 1 so in the many-flavor limit the self-interacting energy error may be neglected. To ensure the B-spline 
grid was sufficiently flne, we compared the trial wave function kinetic and external potential energy before and after 
conversion the B-spline basis set. We also checked the choice of DMC time-step was sufficiently small, the number of 
configurations was suitably large, and the simulation cell size was adequately large. On changing these variables the 
variation in the ground state energy was AE ^ 0.02E^, sufficiently small to allow us to compare the ground state 
energy as the potential strength and gradient expansion coefficient were varied. 

B. Results 

We analyzed a quantum dot containing a MFEG of 12 flavors and 4 bands (shells), containing a total of 96 electrons. 
This was chosen since it had a full shell so is expected to have a zero spin ground-state [13] that can be analyzed with 
the many-flavor functional, was computationally feasible, and contained enough electrons to be in the LDA regime, 
where the many-flavor functional is expected to apply. 

Two different investigations were carried out to probe effects of changing the density gradient, firstly strength of 
the dot confining potential k was changed, and secondly the gradient expansion coefficient ^ was varied. 

1. Varying the external potential strength k 

At the strong external potential k — 8, corresponding to steep gradients. Fig. 3] shows the DFT energy is overes- 
timated compared with the DMC result, indicating that the gradient approximation is not applicable and that the 
next order term in a gradient expansion is negative. Fig. [5] shows that the DFT density profile underestimates the 
true density towards the center of the dot and overestimates density in the outer regions, indicating that the DFT 
functional does not favor steep enough gradients. This is consistent with the next term in the gradient expansion 
being negative. The breakdown corresponds to a coefficient of a « 1.8 in (jmax = a{haQ ^^^)n^/^, close to the a « 4 
which corresponds to the maximum contribution to the interacting energy. 

At the intermediate potential fc = 1 the DFT and DMC estimates of energy and the density profile agree, in 
this region the gradient approximation applies. The DFT density profile shows a slight over-density at the center. 
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FIG. 4: In the upper plot the crosses and solid line show the difference between DFT (-Bdft) and DMC (E) energies with 
varying external potential strength k, if agreement were exact, points would lie along the horizontal dotted line. In the lower 
plot the circles and solid line show the maximum density gradient of the dots with varying k, and the squares and dashed line 
the gradient at which the theory breaks down. 

consistent with self- interaction energy being included in the DFT calculation. At the weak potential fc = 0.1 electron 
densities are low meaning the homogeneous interacting energy is outside of its region of applicability (see Fig. [2]), 
therefore the DFT energy is an overestimate. 

2. Varying the gradient term coefficient ^ 

Fig. [S] shows results of simulations on dots, chosen to have a potential strength k — 1, which is at the center 
of agreement of the previous results. The best agreement between the DFT and DMC ground state energy is at 
^ ~ 0.9. This is in good agreement with the expected ^ = 1, the difference may be due to systematic errors such as 
the self-interacting energy or higher order gradient terms. As expected, the energy is overestimated for dots with too 
large a gradient expansion term, and underestimated for dots with too small a gradient correction term. 

The maximum gradient seen in the dot density profile decreases as ^ increases (see Fig. [5]). The dot becomes more 
spread out so the external energy Ecxt increases whilst the total electron-electron Coulomb energy E^-c decreases. 
Overall the total DFT energy increases. Three quantum dot electron density profiles for gradient term coefficients 
^ = 1,2, and 3 are shown in Fig. [71 Compared with the dot calculated with ^ = 1, the dot generated with no energy 
penalty for gradients, ^ = 0, has a high central and low outer density showing that it has a higher gradient in the 
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FIG. 5: The lower panel shows the density profile of quantum dots estimated using both DFT and DMC at external potential 
strengths of fc = 1 and fc = 8. The difference between the DFT and DMC results at fc = 1 and fc = 8 is shown in the upper 
panel. The DMC statistical error is less than the size of the points. 



density. Conversely dot with increased energy cost for gradients, ^ ~ 2, has a more shallow profile. 

The density profiles seen in Fig . [5] and Fig. [7] can be further analyzed in light of other theoretical studies of 
quantum dots reviewed in Ref. |47| . The density profile calculated using the many-flavor functional is not flat at 
the center, but instead has correlation-induced density inliomogeneity evidenced by a characteristic minimum in the 
density at r « 2aQ. The intermediate density regime in which this occurs is consistent with the strong correlations 
causing a minimum in the total many-flavor energy density 0]. It is also akin to the intermediate density regime 
seen in other quantum dot systems [47l. |5 3 . Iot . [741 . in the high density limit the quantum dot has properties like a 
Fermi liquid with de- localized electrons [43, [S^, uM i whereas in the low density limit the electrons become crystalline 
[47I . [H, [H, [tI, [73, [2I] inside the dot. As the many-flavor functional was successful in predicting correlation-induced 
inhomogeneities, it could be used to investigate other many-flavor quantum dot effects including the Kondo effect in 
multi- valley semiconductors [zl, [s^l , the reduction of valley degeneracy of coupled quantum dots ^ 
harmonically trapped cold atoms with an additional quantum number denoting energy level [ll[, [T^.]l3 



VI. CONCLUSIONS 



We have computationally verified the theory of the MFEG presented in 0] using QMC simulations. In a homoge- 
neous system, DMC estimates for the ground state energy are consistent with theory and the theoretically estimated 




FIG. 6: The upper panel shows, for dots with external potential strength k = 1, the density profile maximum gradient as a 
fimction of ^. The central panel shows the variation of external potential DFT energy based on the primary j/-axis using square 
points and the solid line, the secondary y-axis shows electron-electron DFT energy using crosses and the dashed line. The lower 
panel solid line shows the variation of DFT ground state energy with ^, and the horizontal dotted line the ground-state energy 
predicted using DMC from the ^ = 1 trial wave function. The DMC statistical error is less than the size of the points. 



density range over which the theory apphes is consistent with numerical results. The applicable density for Bi2Te3 
{i^ = 6) corresponds to a charge carrier density between 7 x 10^^ cm^"^ and 2 x lO^*' cni^^. 

The density response function for a MFEG with 24 flavors was found using three methods: density modulation 
predicted by VMC, and the variation in ground state energy predicted by VMC and also by DMC. The two VMC 
results underestimated the response 1/e, but the DMC results agreed with theory and could distinguish between the 
exact and a useful approximate expression for polarizability. 

We used a many-flavor functional including a local gradient approximation in DFT calculations of large quantum 
dots. The DFT calculation estimated the ground-state energy and wave function, which were verified by a DMC 
calculation. We found the high gradient breakdown of the expansion was at ^max ~ l.8{haQ ^^^)n^/^, the low 
gradient breakdown was consistent with the homogeneous MFEG lowest applicable density, and that the gradient 
expansion was applicable in the intermediate regime. The many-flavor functional, used as part of DFT calculations, 
could be a useful tool for analyzing other multi-valley semiconductor systems. 
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density profiles from the ^ = 1 DFT density profile. The lower panel shows DFT density profiles for ^ = 0, ^ = 1 (solid line, 
triangles), and ^ = 2inaA; = l dot. 
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